Density matrix numerical renormalization group for non-Abelian symmetries 
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We generalize the spectral sum rule preserving density matrix numerical renormalization group 
(DM-NRG) method in such a way that it can make use of an arbitrary number of not necessarily 
Abelian, local symmetries present in the quantum impurity system. We illustrate the benefits of 
using non-Abelian symmetries by the example of calculations for the T-matrix of the two-channel 
Kondo model in the presence of magnetic field, for which conventional NRG methods produce large 
errors and/or take a long run-time. 
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I. INTRODUCTION 

Quantum impurity models play a crucial role in our un- 
derstanding of strongly correlated systems: They appear 
in the description of correlated mesoscopic structures^ 
they show up in molecular electronics, and many of the 
properties of correlated bulk systems can also be ac- 
counted for using self-consistent quantum impurity mod- 
els within the dynamical mean field approach.— Despite 
the lot of interest and the large amount of effort in- 
vested in understanding these models, we have, unfor- 
tunately, very limited tools to describe quantitatively 
the general properties of a generic quantum impurity 
model. For some of the quantum impurity models Bcthe 
AnsatZ)^ conformal field theory^ bosonization^ per- 
turbative calculations,—^ or a Fermi liquid theory^ can 
provide a satisfactory explanation. However, these re- 
sults are usually restricted to some regions in the pa- 
rameter space. Therefore, even today, the most reliable 
method to obtain accurate information on a generic quan- 
tum impurity model over the whole parameter space and 
for any frequency is Wilson's numerical renormalization 
group method (NRG) 1 ^ originally developed for the one- 
channel Kondo model (1CKM). 

Apart from being extended to compute dynamical 
properties ; 11 ' 12 ' 13 ' 14 Wilson's method has been used in 
its original form for a longtime, and it is only recently 
that this method has been further developed. First, Hof- 
stetter realized that a density matrix NRG (DM-NRG) 
procedure needs to be introduced in certain cases to avoid 
spurious results.— The method of Hofstetter, however, 
applied the original truncation scheme of Wilson, and 
therefore has not conserved spectral weights. A remark- 
able development of the NRG scheme has been the intro- 
duction of a complete basis set by Anders et al. ^ which 
was used to develop the time-dependent DM-NRG al- 
gorithm, and it also led to the development of spectral 
sum-conserving DM-NRG algorithms . 17 ' 18 In the work of 
Weichselbaum and von Delft, the spectral sum conserv- 
ing method has been hybridized with a matrix product 



state approach.— 

Using symmetries is an important element of NRG. 
For the simplest models, it is usually sufficient to use 
Abelian symmetries. However, for the more interest- 
ing two- and multi-channel models it is crucial to ex- 
ploit symmetries, since the computational effort needed 
increases rapidly with the number of electron channels, 
and it is very important to keep a sufficiently large num- 
ber of levels in the DM-NRG procedure to achieve good 
accuracy. Despite the amazing development discussed in 
the previous paragraph, a general framework in which 
the advantages of non-Abelian symmetries are exploited 
in conjunction with DM-NRG was still missing, apart 
from the generalization of DM-NRG for the case of only 
one SU(2) symmetry^ The aim of the present paper is 
to develop such a scheme and to show, how the spectral 
sum-conserving DM-NRG method can be used in combi- 
nation with non-Abelian symmetries. For this purpose, 
we derive a very general recursion formula for the reduced 
density matrix. We shall show how this formula can be 
used to evaluate the spectral function of any retarded 
Green's function of the form 



G£ flt (t) = -i<[A(t),iJ(o)fc>e(f). 



Here A and B are any kind of local fermionic or bosonic 
operators acting at the impurity site, and [A, B]^ = 
AB — £BA denotes the commutator (£ = 1) and the 
anticommutator (£ = — 1) for bosonic or fermionic oper- 
ators, respectively. In Eq. |T]), (...) denotes the average 
with the equilibrium density matrix, (...) = Tr{. . . g}. 
Expressions for static quantities shall also be derived. 

Due to the general recursion relation mentioned above, 
we were able to perform a DM-NRG calculation indepen- 
dently of the type and number of symmetries considered, 
and to build a completely flexible DM-NRG code that 
handles symmetries dynamically and "blindly" 2^ 

To demonstrate the benefits of our procedure, we shall 
present numerical results for the two-channel Kondo 
model, which is the most basic non-Fermi liquid quan- 
tum impurity model^ and provides an excellent testing 
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ground for multiple non-Abclian symmetries. Further- 
more, its study is also motivated by its recent mesoscopic 
realization in double-dot systems by Potok et al— As we 
shall see, conventional NRG techniques and DM-NRG 
with Abelian symmetries give rather poor results com- 
pared to a computation performed using the non-trivial 
symmetry structure of the model. 

The paper is organized as follows. In Section |TT] we 
present how one can exploit the internal symmetries of 
the quantum impurity system in the NRG process. In 
Section IlIII the general procedure, which permits the use 
of an arbitrary number of local symmetries, is extended 
to the spectral sum rule preserving DM-NRG algorithm. 
In Section IIVI we provide the general formulas for com- 
puting Green's functions in the DM-NRG framework us- 
ing symmetries. In Section [V] we present our numer- 
ical results for the local fermions and local composite 
fermions spectral functions and for the on-shell T-matrix 
in the two-channel Kondo model which strongly support 
the use of DM-NRG (with the largest possible symmetry 
of the system) as opposed to the use of NRG. In Section 
I VII we sum up what has been presented. In Appendix 
IVIII we give an alternative derivation of the result that 
the reduced density matrix retains its diagonal form in 
the course of the reduction for the special case when the 
symmetry of the system contains only SU(2) groups from 
among non- Abelian groups. 

II. THE ROLE OF SYMMETRIES IN THE NRG 
PROCEDURE 

In the present section we shall briefly discuss, how 
non- Abelian symmetries appear in the NRG calculations. 
However, before setting up the general formalism it is 
useful to discuss the role of non- Abelian symmetries in 
the specific example of the one-channel Kondo (1CK) 
model and consider the general case and the general re- 
cursion relations only afterwards. 

A. Symmetries of the one-channel Kondo model 

In the NRG procedure, Wilson used the following ap- 
proximation for the 1CK Hamiltonian^ 

HlCK = ^ 3 E ^fo,^<7<7'fo,a' 

tr,<7'e{T,J.} 

oo 

+ E E *» (4w + h -°) ■ ( 2 ) 

In this Hamiltonian /J a creates a conduction electron 
of spin a at the impurity site, and thus the first term 
describes the interaction between the impurity spin, S, 
and this local fermion through an exchange coupling, J . 
The dynamics of this local fermion is described by a semi- 
infinite chain, the Wilson chain. Electrons (fermions) 



move along this semi-infinite chain with an exponentially 
decreasing hopping amplitude, t n oc A~™/ 2 , with A the 
discretization parameter.— 

The above model has an SUs(2) symmetry, corre- 
sponding to spin rotations, i.e. it is invariant under uni- 
tary transformations 

H\ck = U s H 1C k U\ , (3) 

where the unitary operator IA S is generated by the total 
spin operators, 

U s = U s {u s ) = e iS ' §T , (4) 

1 OO 

St = S+-J2 E fla^'fn,,'- (5) 

The Hamiltonian Hick is also invariant under the ac- 
tion of SUc(2) rotations in charge space, U c = e %u) " , 
generated by the operators C x = (C+ + C")/2, C y = 
(C+-C-)/2i and C z with 

OO 

1 oo 

C " = 2^ £ (fU^-l), (6) 
n=0/*={U} 

c- = (c+y . 

Since the spin symmetry generators commute with the 
charge symmetry generators, Hick possesses a symme- 
try, SUs(2) x SUc(2)r^ and consequently, the eigen- 
states of the Hamiltonian form degenerate multiplets, 
\i,Q.,Q z ), that are classified by their label i, the spin 
and charge quantum numbers: Q. = {Si,Ci\^ and the 
z-components of the charge and spin operators, Q z = 
{Sf,Cf}, 

In the following, we shall refer to the quantum numbers 
Q as representation indices, while Q z are referred to as 
labels of the internal states of a given multiplct. Note that 
the representation index Q defines the dimension of the 
given irreducible subspace which is dim(z) = dim(Q.) = 
(2 Si + l)(2Cj + 1) in the example above. 

Similar to irreducible subspaces of the Hilbert space, 
operators can be organized into irreducible tensor opera- 
tor multiplets.— One of the simplest examples is provided 
by the impurity spin, from the components of which we 
can form the operator triplet as 

{A m }^|--i=S + ,S^,-i=S-| . (7) 

The components of this triplet transform under spin ro- 
tations as the eigenstates \m) of S z within a spin S = 1 
multiplet, while they are invariant under charge rota- 
tions. This means that A has quantum numbers, Sa = 1 
and C a = 0, and the components of this operator mul- 
tiplct arc labeled by S A = m,(m = 0, ±1), while the 
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internal charge label is trivially C Z A = 0. Similar to the 
eigenstates of the Hamiltonian, in our simple example 
the quantum numbers of an irreducible tensor operator 
B can be organized into a representation index vector, 
b = (Sb 7 Cb) and its components are labeled by b z = 



C%) taking the values S% = -S B , S B + 1,...,S B 
and C B = —C B , —Cb + 1, . . • , C B . A further example for 
an S = 1/2 and C = 1/2 operator is formed by the four 
operators {/J T , -fo,i, /o,t}- 



The Wigner-Eckart theorem^ then tells us that, apart from trivial group theoretical factors (Clebsches), the matrix 
elements of the members of a given operator multiplet and states within two multiplcts, i and j are simply related by 



B 



b,b. 



3, 



—3—3 



B II j) (Q..Q 



(8) 



where (i \\ B \\ j) denotes the reduced (invariant) matrix element of B, and the generalized Clebsch-Gordan coeffi- 
cients are simply defined as 



b,b z ;Q.Q* 



(SiS z \ S B , S B \ SjSj') (CiC*\ C B , C B ; CjC*} 



(9) 



with the usual SU(2) Clebsches^ on the right hand side. This relation is used extensively in the NRG calculations. 

I 



An important property of the unitary transformations 
above is that they are local in the sense that they decom- 
pose into unitary operators which commute with each 
other and act independently at different sites, 



U = U s x U c , 



U, 



(10) 

(11) 

(12) 



This decomposition property is crucial for using symme- 
tries in the NRG calculations. 



B. General local symmetries on the Wilson chain 

Under rather general conditions, the symmetry consid- 
erations above and the NRG procedure can be extended 
to Hamiltonians of the form 



This recursion is depicted in Fig. [T] 

Let us now assume that H as well as every H n is in- 
variant under the group G, i.e. 



U(g)H n U- 1 (g) = H n 



n = 0,1,2,. 



(16) 



holds for every g £ G, with U(g) the appropriate unitary 
operator. Furthermore, let us suppose that G and corre- 
spondingly U can be decomposed into a direct product of 
r subgroups Qj (7 = 1, T), each acting independently 
on every lattice site, 

g = g 1 x g 2 x ■ ■ • x g v , (17) 

r r 

u{ g ) = 11^(57) =nrK,»(<?7), (is) 

7—1 7—1 n 

Note that in the considerations above the subgroups can 
be also finite, and G 1 can represent a crystal field sym- 
metry as well as e.g. the SU(3) group. However, some of 



H = Hq + ^2 ( T n,n+1 + 7~i n+ i) 



(13) 



Here /Hq contains the interaction between the impurity 
and the fcrmionic bath (the site of the impurity is labeled 
by 0), and nearest-neighbors on the Wilson chain are 
coupled through the hopping terms, r nj „+i. The n th on- 
site term H n +i describes local correlations/interactions, 
and it can also account for the absence of the electron- 
hole symmetry. 

The usual NRG solves the model Eq. (fT3"|) by an itera- 
tive diagonalization process. The iteration steps consist 
of diagonalizing the set of Hamiltonians introduced re- 
cursively by 

H = Ho , (14) 

H n+ i = H n + T n ,n+l + Hn+1 ■ (15) 




Tn-i,n 



FIG. 1: Construction of the Hamiltonian on the Wilson chain 
of length N. The impurity sits at the square shaped site 
labeled by 0, dots represent the further sites. At the zero" 1 
NRG iteration Ho is identical with Ho- As the n th site is 
added to the chain, H n is constructed recursively from the 
hopping term (r n ,»+i), from the on-site energy (Tin), and 
from H n -i. 
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the considerations presented in this paper may not apply 
for non-compact groups. 

The above decomposition is not necessarily unique. 
Nevertheless, having obtained a specific decomposition 
of the symmetry, the argument of the previous subsec- 
tion can be repeated with the only difference that now 
r number of quantum numbers classify the irreducible 
subspaccs (multiplcts) of the Hamiltonians H n , 



Q = {Q 1 ,Q 2 ,...,Q r } 



(19) 



and states within the multiplet are then labeled by the 
internal quantum numbers 



Q z = {Q hz ,Q 2 



(20) 



Similar to the simple 1CK example, the dimension of a 
subspace i depends uniquely on its quantum numbers 
i.e. dim(i) = dim(Q ). 

Operators can also be arranged into irreducible tensor 
operators, and an irreducible tensor operator multiplet 
A is correspondingly described by quantum numbers a, 
while members of the multiplet are labeled by a z with a 
and a z being T-component vectors. The Wigner-Eckart 
theorem, Eq. ([5]), carries over to the general case too with 
the slight modification that the Clcbsch-Gordan coeffi- 
cients are now defined as 



q, q: 



a a 



NRG. In his original work, Wilson constructed approxi- 
mate eigenstates of H n for a chain of length n, iteratively. 
However, in each iteration the dimension of the Hilbcrt 
space increases by a factor d, with d the dimension of 
the local Hilbert space at a single site of the chain with 
site label n > 0. Therefore the size of the Hilbert space 
increases exponentially with n, and after a few iterations 
one must truncate it: Some of the states i arc therefore 
discarded (i G D), while other states are kept (i G K), 
and are used to construct approximate eigenstates for 
H n +i- Symmetries are of great value in this diagonal- 
ization procedure: In their presence the H^s are block 
diagonal in the representation indices, and the eigenvalue 
problem can be solved much more efficiently. For some of 
the physical quantities, it is crucial to increase the num- 
ber of kept states as much as possible to achieve good 
numerical accuracy. 

In Wilson's original formulation of NRG, symmetries 
are used in the following way: As discussed above, in 
the n th iteration the eigenstates (multiplets) of H n are 
constructed from the kept multiplcts labeled by u with 
u G K of the (n — l) th iteration, which are approximate 
low-energy eigenstates of H n —i, and from a complete set 
of local states (multiplcts), labeled by fj., that live at the 
n th site. In the following, we refer to these new approx- 
imate eigenstates i as new states, while by analogy with 
DMRG, we shall call the kept states block states or old 
states. 



n (Qi ®r 



C. The role of symmetries in the diagonalization 
procedure 

Before discussing NRG with a complete basis set^ let 
us shortly look into how symmetries are used in Wilson's 



In the presence of symmetries, each new multiplet car- 
ries representation indices Q~ = {Q~}, and states within 
this multiplet are labeled by the internal quantum num- 
bers — {Qj' Z }- Similarly, local states have quantum 
numbers q = {q2} and are labeled by q z = {qZ' z }- To 
construct the approximate eigenstates of H n , we first use 
the Clebsches to build new states from the block (old) 
and local states that transform as irreducible multiplets 
under the symmetry transformations, U(g), 



< -OA/ 



E 



q q z ;Q Q" 



u,Q u ,Ql;fi,q^ql) n _ 



{u G K) 



(21) 



These states shall be referred to as the canonical basis 
from now on. In this basis H n has a block-diagonal struc- 
ture 



i, Qj, Qf | H n | j, Q v Q| 



(* II H n || ]) n Sq^qJq^q. , 



and is diagonalized by a unitary transformation, D\ . 

=£l°!5 '•^X/V'^ .:22; 

i,Q,,Q z ) , (23) 



Ef 



with E™ the eigenenergies of H n . Here O is a block- 
diagonal matrix, and its columns in a given symmetry 
sector are just the eigenvectors of the corresponding sub- 



matrix of (i \\ H n \\ j) in the canonical basis. In the up- 
coming iteration, some of these multiplets shall be kept, 
and form the block states for (n + l) th iteration, while 
others are again discarded. 

In the iteration step outlined above we need the matrix 
elements n (i || H n \\ j) n , with H n = H n -i + T n -i, n + 
Ti. n - The matrix elements of H n -\ are simply 



(i \\H n . x || j) n = K 



1 *u 



Here cj^ c ^l denotes some creation operator multiplet at 

site n — 1 that has quantum numbers c, and ha™ 1 ' 's arc 
the hopping amplitudes between sites n — 1 and n. The 
index a in the equation above labels various " hopping 
operators". To give a simple example, if we treat the 
1CKM using only U(l) symmetries then one has two hop- 
ping operators, a £ {1,2}, corresponding to c["' = * 
(24) and CJ,™ 1 = A , • However, if we use the spin SU(2) sym- 



where u is the state from which state i has been con- 
structed. Similarly, matrix elements of Ti n are given by 



(i II H n 



(25) 



where e™ is just the expectation value of Tt n with local 
states within the multiplet (i. Finally, to compute the 
matrix elements of the hopping t„_i.„, we use the fact 
that by symmetry, T n -i, n can always be decomposed as 



mctry, then a 



land ^ = {/t T)/ t ;T }. 



7~n— l,n 



(26) 



For the reduced matrix elements of r n _i iTl one obtains 
using Eq. (|2"Tj) and the decomposition, Eq. ([26|) . the fol- 
lowing formula, 



Q II II J> „ = E n-l ( U II II ») n-1 (" I" C ^ I" X 

a 

D i K ^c,Q, : Q,,Q^ u ,Q^,q^q^5Q j> Q_5^Qz . (27) 
I 



Here the state i has been constructed from the state u of 
the previous iteration and from the local state [i, while 
j has been constructed from v and v. The functions 



D yjtjCjQ^QpQ^Q^q^q^j denote group theoretical 
factors, 



D (a.ft&.fij.S.'S.' V2j = sg»(C, /*) E E E (QlSS 

C' Qf, , Q* il , g* 



q q~;ytj~) X 



Q Q 2 



cc z ;QQl)(qQ 



V — —V / 



: W (28) 



where the sign function sgn(C, /i) = ±1 arises as one 
commutes the creation operators implicitly present in the 
local state \i over the operator C^ n ~ l \ and it is negative 

r n i 

if C a is a fermionic operator and the local state \i is 
constructed from an odd number of fermions, otherwise 
it is positive. The local matrix elements, (// || cl™' || v) 



are the same for all sites and can easily be determined, 
while || C'™ -1 ! || v) _ 1 can be computed from the 

previous iteration by recursion. In fact, for any operator 
A, acting on sites m < n, and whose matrix elements 
are known in the iteration n — 1, we have the following 
recursion relation 



(i || A || j) n = n _ x < u || ^ || v) n _ x F (a,Q- v Q v Q u ,Q v ,q^ 6 q _^ , (29) 
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where the factor F ^a, Q^, Q~., Q u , Q^, q^j is given by the following expression 

*(«.a.$.ft,a..O ") dim( ^, m(j) £(&&■ 



aa z \Q~.Q z ) x 



E E(<M 



9 <Z Z ;Q Q Z )(Q,Q? 



? q z ;Q Q z ) (Q Q z 



aa z ;QQl) (30) 



One drawback of the algorithm above is that the eigcn- 
states of H n constructed this way do not form a complete 
basis set on the Wilson chain of length n, since states 
descendant from the discarded states of the previous it- 
eration are missing. However, as it was recently shown,— 
one can construct a complete basis of the Wilson chain 
in a slightly different way: Let us consider a chain of 
length N, and construct approximate eigenstates of H n 
that, however, live at all sites of this chain, 



■ Q.tr 



,Q Jt Q* 



(31) 



In this equation e just labels the d N ~ n independent 'envi- 
ronment' states living at the last N — n sites of the chain. 
The internal structure of these environment states is not 
important, only their degeneracy shall play some role. 
The previous iterative construction carries over to these 
states, too. By construction, discarded states (together 
with their environment state) form a complete basis set: 



N 

1 =EEEE 

n=0 ieB e Qz 



(32) 



where i € D refers to the fact that only discarded states 
appear in the sum. In Eq. (|32|) all states are consid- 
ered discarded in the last iteration, n = N. We remark 
that, in reality, the summation in the expression above 
starts only at a value n = M, where the first truncation 
is carried out. Fig. [1 illustrates the structure of this com- 
plete basis. In the formulation of the DM-NRG algorithm 
we shall use several times the completeness relation, Eq. 
(El). 



III. CONSTRUCTION OF THE REDUCED 
DENSITY MATRIX USING SYMMETRIES 

In the DM-NRG procedure on a Wilson chain of length 
N, the equilibrium density matrix is approximated by 



/v 



E e w 



= E 



-0E7 



(33) 



(34) 



Kept states 



Discarded states 




u 

a 



< 



M 



M+l 



M+2 



FIG. 2: (color online) A complete basis of a Wilson chain rep- 
resented as the exponentially increasing number of energy lev- 
els belonging to the successive iterations. Continuous/dashed 
lines represent kept, low-energy/discarded, high-energy levels, 
respectively. For the consecutive iteration steps the distances 
between the levels illustrates how the energy resolution of 
NRG gets exponentially refined. 



with (3 = 1/ksT the Boltzmann factor and 

JV 

n=0 Q*,i 



(35) 



the partition function. In Eq. (f3"5|) the factor d ac- 
counts for the degeneracy of the environment states in 
iteration n, i.e. for the local degrees of freedom at sites 
m > n. Since the eigenenergies do not depend on the in- 
ternal quantum numbers, the expression for the partition 
function can be simplified to 

N 

Z = J2Y, dim W e~ 0E? d N - n , (36) 

(37) 



n— i 



dim(<) = [] dim(Q7) , 

7 = 1 



with dim(Qj) being the dimension of an irreducible sub- 
space having the representation index Qj. 

The concept of the reduced density matrix^ arises nat- 
urally as one starts to calculate Green's functions with 
NRG. More precisely, the quantity that shows up in the 
calculations is the truncated reduced density matrix, de- 
fined as 



R [n\ = Tr 
{en} 



E 



(38) 
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where one traces over environment states e n at sites m > 
n. This truncated reduced density matrix clearly satisfies 
the recursion relations 



R [N] = g [N] 



R l 



Tr 

site n 



-/3EV-' 1 jN-n+1 



(39) 



Note that the environment variable is missing in the first 
term of the second expression since it has been traced 
over. The first term accounts for the contribution of dis- 
carded states, while the second term has matrix elements 
between the kept states only. 

To construct the matrix elements of R^ we first show 
by induction that R^ is scalar under symmetry opera- 
tions. This is clearly true for the first term in Eq. ([59")) . 
To show that the second term is also invariant, we simply 
need to use the locality property of the symmetry trans- 
formations, i.e. that on the first n sites U(g) = L(g)V(g), 
where L{g) transforms the local states on site n, while 
V(g) transforms states at sites < m < n— 1, and L and 
V commute with each other. Therefore, because of the 
trace of operators being invariant under cyclic permuta- 
tions, we have 



{LVR [n] V+L+] = Tr {vR [n] V+\ 

I J site n I J 

= V Tr \R [n] \v + = U Tr \ R [n] \ U + . (40) 

site n I J site n I J 



However, since U B> n ' U + = Ry 1 ' by assumption, we im- 
mediately obtain that 



Tr 

site 



fijW) = U Tr IrW)u+ , (41) 

7i I J site n I J 



implying that 



UR [n ~ 1] U+ = R [n ~ 1] 



(42) 



for : 1 , too. This equation means that is scalar 
and therefore, by the Wigner-Eckart theorem we have 



i? 1 



j,Q„Q' 



Ji \\R [n] II 3)Jq.,qJqw (43) 



The matrix elements of (i \\ R^ |j j) between discarded 
states simply derive from the first term in Eq. ([59"]) . To 



perform the trace in Eq. (|39)) and to construct the ex- 
plicit relation between the kept matrix elements of R^ n ~^ 
and some more work is needed. First, we rotate 
to the canonical basis 

Ji || RW || ~j) n = 0Wji || RW II j)n (O-^ .(44) 
Then, using the fact that Tr is diagonal in the 

site n 

symmetry quantum numbers and labels, we can trace 
over the local states at site n using the recursion relation 
Eq. (|2"Tj) to obtain the matrix elements between the kept 
states u,v £ K as 



i || i?^ 



dim I 



i,3,q..,v 



dim (u) n 



(i\\R ln] \\j)jQ„Q,- (45) 



Here the tilde over the sum indicates that in the sum- 
mation over i and j only those states are considered which 
have been constructed from u (i <= u) and v (j <= v) in 
Eq. (|2ip , respectively. This is a very powerful expression, 
which applies to essentially any type of symmetry. 



IV. SPECTRAL FUNCTION COMPUTATION 

The quantity that describes the linear response of 
a static system to a time-dependent perturbation and 
which is to be calculated in the DM-NRG framework is 
the retarded Green's function. For two irreducible tensor 
operators it is defined as: 



G 1 



(t) 



Aa,a(t),Bl b (0) 



0(t). (46) 



By symmetry, however, this Green's function is non-zero 
only if the spectral operators A and B transform accord- 
ingly to the same representation i.e. 



(*) 



Bt (*) $a,b $a z ,b z 



(47) 



and it is independent of the value of a z = b z . Note 
that in this expression the representation index b and its 
labels b z are the quantum numbers that characterize the 
operator B and not B* . 

In the reduced density matrix formalism we can gener- 
alize the procedure outlined in Ref. [l6j even in the pres- 
ence of non-Abclian symmetries to obtain the following 
form for the Laplace transform of the Green's function 
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^A,Bf( z ) 



JV 

E 

n=0 i€DJ< (j,fc)£(K,K) 



E E 



(fc || At || j) * n <fc || St || i) n dim(fc) t n ( j || g || k)nn{i || ^ || fc);dim(») " 

z- \{Ef + EV-) -El dim(o) 



dim(a) 



(48) 



Remarkably, this formula contains exclusively the re- 
duced matrix elements and the dimensions of the various 
multiplets. Here the second sum is over all the multiplets 
k of the given iteration subject to the restriction that 
j, k do not belong to kept states at the same time and no 
summation is needed for states within the multiplets. In 
Eq. (|4*5j) dim(a) = Jl^i dim(a 7 ) is the dimension of the 
operator multiplet Ag,,a ■ We note that the irreducible 
matrix elements of are identical with the original 
ones since Ry 1 ' is invariant under all symmetry trans- 
formations, i.e. it is a rank object with respect to all 
symmetries. Eq. (|48|) explicitely shows that unless A and 
B have the same quantum numbers (a=b), Bt (z) = 0. 

V. NUMERICAL RESULTS 

In this section we show the advantages of using DM- 
NRG as opposed to NRG, and illustrate the benefits of 
using non-Abelian symmetries by applying DM-NRG to 
the two-channel Kondo (2CK) model. This model is ex- 
citing in itself as it possesses a non-Fermi liquid type of 
fixed point and it provides the simplest descriptions of 
the double dot system used recently to realize the 2CK 
stated This 2CK state is very fragile and in a magnetic 
field the difference between the NRG and DM-NRG re- 
sults are substantial. 

In the Wilson approach, the 2CKM is described by the 
following Hamiltonian, 

H2CK = 2 S E ^ a E fo,u,<r a aa' fo,a,a' + 
QS{1,2} <t,<t'G{T,1} 

OO 

+ E E E *« (fL,Jn,a,„ + h.C.) , (49) 
n=0 qS{1,2} <re{T,i} 

where we have introduced a G {1,2} for labeling the 
two types of electrons. This additional channel label is 
the only difference compared to the one-channel Kondo 
Hamiltonian (cf. Eq. ©). 

In this model, the number of carriers is conserved in 
both channels corresponding to a Uci(l) x Uc2(l) sym- 
metry. However, due to the presence of electron-hole 
symmetry, these charge symmetries are augmented to 
SU(2) symmetries, and the Hamiltonian above is also 
invariant under SUg(2) x SUci(2) x SUc*2(2) transfor- 
mations. On the other hand, one can solve this Hamilto- 



nian using exclusively U(l) symmetries, Us(l) xUci(l) x 
Uc2(l)- This model is thus ideal for testing our flexible 
methods. 

If a local magnetic field is coupled to the impurity spin 
through a term gfj,BBS z , then from among the total spin 
generators (Eq. |(SJ)) solely will commute with the 
Hamiltonian. That is, the spin SU(2) symmetry of the 
system reduces to U(l). Therefore, in a magnetic field we 
can either use the symmetry Us(l) x SUci(2) x SUc2(2) 
for our calculations, or restrict ourself to U(l) symme- 
tries only: U s (l) x U C i(2) x U C2 (2). 

As a test, we computed the retarded Green's function, 
G 1 } t (uj) and the corresponding spectral function 

*0,c«,T'-'0,q,t 

Pf 0aT (uj) both in the presence and in the absence of 
magnetic field. All numerical results presented were ob- 
tained at zero temperature, and the dimensionless cou- 
plings were J7q = 0.2 for both channels and all runs. The 
discretization parameter A = 2 was used in all cases, 
and for each symmetry combinations we have retained a 
maximum number of 1350 multiplets in each iteration. 

In Fig. [3] we show data for the local fermion's spec- 
tral function in the absence of magnetic field obtained 
through the NRG and the DM-NRG approaches using 
the two symmetry groups mentioned above. The Kondo 
scale Tk in Fig. [3] is the scale at which the 2CK state 
forms, and it is defined as the frequency where the T- 
matrix of the 2CK model drops to half of its value as- 
sumed at to = at the 2CK fixed points 

The first important test is the fulfillment of the spec- 
tral sum rules. These are always satisfied in the DM- 
NRG calculations independently of the symmetry group 
used, whereas the NRG data violate the sum rule to 
over 15% if the number of kept multiplets is 1350 cor- 
responding to m 7 x 10 4 states. Fig. [31(6) shows that 
the expected t/uj behavior around the 2CK fixed point 
is nicely recovered by both methods, but a sufficiently 
large number of multiplets must be kept in the DM-NRG 
approach, meaning that in this case the larger symme- 
try group must be used. Fig. [31(c) demonstrates that, 
in spite of fulfilling the spectral sum rules still the DM- 
NRG data do not show the expected asymptotics for low 
frequencies if the number of multiplets kept is not suf- 
ficient. It is quite remarkable that only the DM-NRG 
procedure using non-Abelian symmetries was able to get 
close to the exact value of the spectral function at w = 0, 
P/o, Q . T (^-0) = 0.25. 
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Spectral weights 



B/T K =0 




FIG. 3: (color online) Dimensionless spectral function of /o,i,t 
normalized by D, the bandwidth cut-off, as a function of 
uj/Tk in the absence of magnetic field obtained with DM- 
NRG and with NRG using the symmetries: Us(l) x SUci (2) x 
SU C2 (2) and Us(l) x Uci(l) x U C2 (1). (a) Comparison be- 
tween the spectral weights of the DM-NRG and NRG results: 
DM-NRG fulfills the sum rule entirely even when the used 
symmetry group and therefore the number of kept states is 
largely reduced. NRG violates the sum rule to over 15% if the 
number of kept states is ~ 7 x 10 4 in each iteration, (ft) The 
same spectral functions as a function of \Juj/Tk- If a suffi- 
cient number of states is kept, i.e. when using larger symme- 
try groups, the expected ^/uj behavior around the 2CK fixed 
point is nicely recovered, (c) The same spectral functions on 
a logarithmic scale. 



The presence of magnetic field generates a new scales 



Th = C h — 
-Ik 



(50) 



where we have fixed the somewhat arbitrary constant 
to Ch ~ 60.— This scale is usually referred to as the 
renormalized magnetic field acting on the impurity, and 
below this scale the non-Fermi liquid physics is destroyed. 

In the presence of magnetic field the sum rule is vio- 
lated by the NRG approach to a different extent in the 
positive and negative frequency ranges, which leads to 
jumps at u> = in the spectral functions (see Fig. Q|, 
while this problem is absent in the DM-NRG approach. 

Spectral functions display universal scaling in the 
vicinity of Th<^ In Fig. [5] we show how the spectral func- 
tions of the composite fermion operator, F± = /J 
can be scaled on top of each other using the scale Th- Al- 
though, this collapse can be obtained in both approaches, 
there is an sw 20% jump at oj = in the NRG results 
while the DM-NRG results are continuous there. It is 
possible to eliminate the jump in the NRG results by 
determining the phase shifts from the energy spectrum 



3 0.4 




NRG, U s (1) x SU C1 (2) x SU C2 (2) 

— DM-NRG, U s (1)xU c ,(1)xU c2 (1) 
.- DM-NRG, U s (1)xSU cl (2)xSU C2 (2) 



Spectral weights 

DM-NRG: 1.00000 
NRG: .80775 



(a) 



60 



FIG. 4: (color online) Dimensionless spectral function of /o,i,f 
normalized by D, the bandwidth cut-off, as a function of 
uj/Tk in the presence of magnetic field obtained with DM- 
NRG and with NRG using the symmetries: U s (l) x SU C i (2) x 
SU C2 (2) and U s (l) x U C i(l) x U C2 (1). (ft) On a smaller scale 
at u) — we show the smoothness of the DM-NRG data using 
both groups and the jump in the NRG results using the larger 
group. 
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FIG. 5: (color online) Dimensionless spectral function of the 
J.- (a), (c) and the |-spin (ft), (d) components of the local com- 
posite fermion operator normalized by D, the bandwidth cut- 
off, for sufficiently small values of B as a function of tu/Th 
scaled on top of each other using NRG and DM-NRG together 
with the group U s (l) x SU C i(2) x SU C2 (2). 



with high precision, but even after these corrections, the 
results continue to violate the sum rule and numerical er- 
rors for low frequencies remain of the same size as before. 
Moreover, the 'universal' curve obtained by conventional 
NRG is clearly incorrect and different from the DM-NRG 
result. Also, if we try to compute the imaginary part of 



the on-shell T-matrix of the 2CK model where the local 
composite fcrmion spectral functions for both f- and |- 
spin components have to be summed up, we end up with 
large numerical errors in the NRG results (see Fig. [5}, 
while DM-NRG provides satisfactory results even in this 
case. It is thus clear from these examples that the DM- 
NRG method together with the use of lots of symme- 
tries produces much more reliable results than NRG with 
non-Abelian symmetries or DM-NRG with only Abclian 
symmetries, and its use is needed to do computations for 
more delicate quantum impurity models. 



VI. CONCLUSIONS 

To summarize, in this paper, we have shown how the 
recently developed spectral sum-conserving DM-NRG 
methods can be used in the presence of any number and 
type of non-Abelian symmetries. The most important re- 
sult of this paper is a very general and simple recursion 
relation for the truncated reduced density matrix, which 
enables one to compute the spectral properties of any lo- 
cal correlation function in the presence of non-Abelian 
symmetries. The expressions derived hold for almost 
any symmetry, including non-Abelian finite groups, point 
groups, SU(N) groups, and, of course, Abelian groups 
such as U(l) or the parity. Using these symmetries, re- 
duces considerably the time needed for the computations 
and enhances significantly the accuracy of the numerical 
results for dynamical quantities. The general formula- 
tion presented in this paper also allowed us to construct 
a flexible DM-NRG code where symmetries are handled 




FIG. 6: (color online) Imaginary part of the on-shell T- 
matrix using DM-NRG and NRG together with the group 
Us(l) x SUci(2) x SUc2(2) for various magnetic field values 
as a function of lo/Tr-- 
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dynamically, and which is able to learn and handle es- 
sentially any type of symmetry^ 

We demonstrated the advantages of the generalized 
method by applying it to the two-channel Kondo model 
in a magnetic field. The presence of magnetic field makes 
this model very challenging from the point of view of 
NRG calculations. Wc have carried out calculations for 
the local fermion's and for the local composite fermion's 
spectral function at zero temperature and in a finite mag- 
netic field. In conventional NRG, there is always a jump 
at uj — > between the positive and negative frequency 
parts of the spectral function. In addition to being spec- 
tral sum-conserving, the use of DM-NRG almost com- 
pletely eliminated this jump and made it possible to com- 
pute the universal cross-over curves in a magnetic field. 
Moreover, the DM-NRG approach used with larger sym- 
metry groups has provided substantially better results for 
the spectral functions: In the limit w^Owe needed to 
use non-Abelian symmetries in the DM-NRG approach to 
recover the expected power law behavior known from con- 
formal field theory. Also, the universal crossover curve 
and the shape of the the peak at the renormalized mag- 
netic field Th oc B 2 /Tk was much better resolved by 
DM-NRG than by NRG, and thus DM-NRG with non- 
Abelian symmetries lead to much more reliable results for 
the T-matrix (see Fig. [5]) than NRG with non-Abelian or 
DM-NRG with Abelian symmetries. We emphasize that, 
to obtain the u> = value of the spectral function cor- 
rectly as well as its proper scaling behavior, we needed 
to use non-Abelian symmetries. 

The method presented here thus opens up the possi- 
bility to carry out very accurate spectral sum-conserving 
DM-NRG calculations for multichannel systems, such 
as multi-dot devices, and to perform reliable DM- 
NRG-DMFT calculations for multichannel lattice mod- 
els. Combined with the matrix product state approach, 
it might also provide a way to use methods applied in 
DMRG to improve the high frequency resolution of NRG. 
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VII. PROOF OF THE DIAGONAL FORM OF 
THE REDUCED DENSITY MATRIX FOR SU(2) 
SYMMETRIES 

In this appendix we prove, in a different way from how 
it was demonstrated in the main part of the article, the 
general theorem for the diagonal form of the reduced den- 
sity matrix presented in Section Hill in the case when only 
SU(2) local symmetries are involved. 

The proof for the diagonal form of the reduced den- 
sity matrix goes by induction for the iteration steps. At 
the last iteration, by construction the reduced density 
matrix B} N ^ is scalar under the local symmetry group 
(cf. Eq. ([39]) ). Due to this invariance by very general 
considerations R} N ^ can be written as a sum over tensor 
products of irreducible tensor operator components, and 
is of the form 



following matrix elements 



^ = E( T « oc] )! t .®( T « JV ~ 11 )tt- 



(51) 



(« II II «> = 



<x,q ,m 



x (u || Tl N -V || V ) (jn || rif°°] || M ) 
x (Q Q z \tf Q Q 



E (iX^zX) ■ ^ 



9f, = — '?,, 



In Eq. (|52|) we have applied the Wigner-Eckart the- 
orem. The sign factor sgn(., .) depends on the number 

of fermionic operators used in the construction of T„ ^ 
and the local states. In Eq. ([52]) in the last sum only the 
terms with t z = give non-vanishing contributions and 
the sum can be reduced to 



where (tJ, ,0uj J ^ and (IT _±J J ( 

are irreducible tensor 

operator components of the same rank t acting on the lo- 
cal vector space at site N and on the rest of the Wilson 
chain, respectively and a labels all possible tensor opera- 
tors. This special form is a consequence of the invariance 
of RV*] under the local symmetry transformations.— 

To obtain R^k^ we nave to trace over the local ba- 
sis states, | \i, q^, ?*^}> that is we have to compute the 



E (qq Z \tOqq z ) = [2q + 1) 5 



t,0 ■ 



(53) 



Eq. (|53|) implies that the only non-vanishing contribu- 
tions are those corresponding to t — t z — 0, i.e. to scalar 



rI' ocl and Tl N ~ 



Therefore the reduced density matrix 



R^- 1 ] is diagonal in the representation indices. The in- 
duction towards smaller iteration steps goes recursively 
the same way. 
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